import scvelo as scv
scv.logging.print_version()
scv.settings.verbosity = 3 # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True # set max width size for presenter view
scv.set_figure_params('scvelo') # for beautified visualization
adata = scv.read('path/file.loom', cache=True).adata = scv.datasets.pancreas()
adata
Names of observations and variables can be accessed via adata.obs_names and adata.var_names, respectively. AnnData objects can be sliced like dataframes, for example, adata_subset = adata[:, list_of_gene_names].
scv.pl.proportions(adata)
フィルタリングと正規化は、スプライスされた/されていないカウントにも同じように適用されます。
これらはすべて単一の関数 scv.pp.filter_and_normalize にまとめられており、以下のコマンドが含まれます。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
#scv.pp.filter_genes(adata, min_shared_counts=20)
#scv.pp.normalize_per_cell(adata)
#scv.pp.filter_genes_dispersion(adata, n_top_genes=2000)
#scv.pp.log1p(adata)
scv.pp.momentsにまとめられています (内部でscv.pp.pcaとscv.pp.neighborsを利用)。scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
逆に、負の速度は遺伝子がダウンレギュレーションされていることを示している。
完全な動的モデルの解は、mode='dynamical'を設定することで得られます。このモードは事前に scv.tl.recover_dynamics(adata) を実行しておく必要があります。このモデルについては次のチュートリアルで詳しく説明します。
scv.tl.velocity(adata)
計算された速度は、カウント行列と同様に adata.layer に保存されます。
approx=Trueを設定することで、縮小されたPCA空間上で計算することもできます)。scv.tl.velocity_graph(adata)
マルコフ遷移行列には scv.utils.get_transition_matrix でアクセスできます。
前述したように、velocityを低次元空間に埋め込むために、scv.tl.velocity_embeddingを使用して得られた遷移確率に対する平均遷移を計算します。
scv.tl.terminal_states で参照できます。最後に、速度は任意の低次元軸に投影され、以下のいずれかの方法で可視化されます。
scv.pl.velocity_embedding,scv.pl.velocity_embedding_grid,scv.pl.velocity_embedding_stream.デモデータには既に計算済みの UMAP 埋め込みと注釈付きクラスタが含まれています。
scv.tl.umap と scv.tl.louvain で取得できます。詳細については、scanpyのチュートリアルを参照してください。basis='umap'とcolor='cluster'を使用するように設定されています。scv.pl.velocity_embedding_stream(adata, basis='umap')
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
これはおそらく最も重要な部分です。つまり、生物学的な結論を投影されたvelocityに限定するのではなく、推定された方向が特定の遺伝子にどのように支持されているかを理解するために、相図 (phase portraits)を介して個々の遺伝子のダイナミクスを調べるべきです。
spliced vs. unspliced phase portrait をどのように解釈するかについては、こちらのgif を参照してください。
さて、いくつかのマーカー遺伝子の相図を scv.pl.velocity(adata, gene_names) または scv.pl.scatter(adata, gene_names) で可視化してみましょう。
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
逆に、負の速度は遺伝子が発現抑制されていることを示している。
例えば、Cpeは発現上昇したNgn3(黄色)からPre-endocrine(オレンジ)からβ細胞(緑)への方向性を説明し、Adkは発現抑制されたDuctal(濃い緑)からNgn3(黄色)から残りの内分泌細胞への方向性を説明している。
scv.pl.scatter(adata, 'Cpe',
color=['clusters', 'velocity'],
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
私たちは、結果として得られるベクターフィールドと推定される系統を説明するのに役立つ遺伝子を特定する体系的な方法を必要としています。そのために、どの遺伝子がクラスター特異的な微分速度発現を持ち、残りの集団と比較して有意に高いか低いかを検定するモジュールを提供します。
scv.tl.rank_velocity_genes モジュールは differential velocity t-test を実行し、各クラスタの遺伝子ランキングを出力します。min_corr)。scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
df.head()
kwargs = dict(frameon=False, size=10, linewidth=1.5,
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
scv.pl.scatter(adata, df['Ngn3 high EP'][:5], ylabel='Ngn3 high EP', **kwargs)
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs)
The genes Ptprs, Pclo, Pam, Abcc8, Gnas, for instance, support the directionality from Ngn3 high EP (yellow) to Pre-endocrine (orange) to Beta (green).
RNA速度によって検出される細胞周期は、細胞周期スコア(位相マーカー遺伝子の平均発現量を標準化したスコア)によって生物学的に肯定されます。
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
周期性のある管状細胞 (Ductal cells)については、SとG2Mの位相マーカーを介してスクリーニングすることができます。上のモジュールはまた、位相マーカー遺伝子のランク付けやソートに使用できるスピアマン相関スコアを計算しており、それらのphase portraitsを表示するために使用することができます。
s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)
s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).index
g2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).index
kwargs = dict(frameon=False, ylabel='cell cycle genes')
scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)
特にHellsとTop2aは、循環する前駆体のベクトル場を説明するのに適している。Top2aは、G2M段階で実際にピークを迎える直前に高い速度が割り当てられる。そこでは、負の速度は、その後のダウンレギュレーションの直後に完全に一致する。
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
さらに2つの有用な統計量があります。
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
これらの結果から、細胞の分化のペースが遅い場合と早い場合、またその方向が不明な場合があることが明らかになった。
クラスターレベルでは、細胞周期の終了後(Ngn3低EP)には分化が実質的に加速され、ベータ細胞の産生中にはペースを維持し、アルファ細胞の産生中にはペースが低下することがわかった。
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
scv.pl.velocity_graph(adata, threshold=.1)
さらに、このグラフは、指定された細胞の子孫/先祖を描くために使用することができます。ここでは、前内分泌細胞は、その潜在的な運命までトレースされています。
x, y = scv.utils.get_cell_transitions(adata, basis='umap', starting_cell=70)
ax = scv.pl.velocity_graph(adata, c='lightgrey', edge_width=.05, show=False)
ax = scv.pl.scatter(adata, x=x, y=y, s=120, c='ascending', cmap='gnuplot', ax=ax)
グラフからroot細胞からの分布を推定した後、root細胞からグラフに沿って歩いた後、その細胞に到達するまでにかかる平均ステップ数を測定することで計算されます。
diffusion pseudotimeとは逆に、root細胞を暗黙的に推定し、similarity-based diffusion kernelの代わりに有向速度グラフに基づいています。
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
scv.tl.paga(adata, groups='clusters')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
これは、左/行から右/列へと読み取ることができ、例えば、DuctalからNgn3 low EPへの信頼性ある遷移を割り当てることができます。
この表は、UMAP上に埋め込まれた有向グラフでまとめることができます。
scv.pl.paga(adata, basis='umap',
size=50,
alpha=.1,
min_edge_width=2,
node_size_scale=1.5)